import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import seaborn as sns
import random
from statsmodels.graphics.tsaplots import plot_pacf
data = pd.read_excel(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\2019_PM25_1g.xlsx", skiprows = 5)
data.rename(columns = {'Kod stanowiska':'Date'}, inplace = True)
# data.set_index('Date', inplace = True)
# data.index = pd.to_datetime(data.index)
first_column = data.iloc[:, 0]
data = data.drop(data.columns[0], axis = 1)
data.columns = data.columns.str.replace('-PM2.5-1g', '')
data.insert(0, 'Kod stanowiska', first_column)
data.rename(columns = {'Kod stanowiska': 'Date'}, inplace = True)
data.head(5)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\4022153991.py:7: FutureWarning: The default value of regex will change from True to False in a future version.
data.columns = data.columns.str.replace('-PM2.5-1g', '')
| Date | DsDusznikMOB | DsJaworMOB | DsJelGorOgin | DsWrocAlWisn | DsWrocWybCon | KpBydPlPozna | KpMogiNowMOB | KpToruDziewu | KpWloclOkrze | ... | SlBielPartyz | SlKatoKossut | SlZlotPotLes | WmElbBazynsk | WmGoldUzdrowMOB | WmOlsPuszkin | WpKaliSawick | ZpSzczAndr01 | ZpSzczBudzWosMOB | ZpSzczPils02 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | 33.40530 | 51.38780 | 118.7730 | 102.0900 | 107.0610 | 64.1177 | NaN | 24.0030 | 51.317 | ... | 110.1990 | 76.7306 | 26.3444 | 34.3706 | 14.9449 | 40.9183 | 75.2000 | NaN | NaN | 73.93500 |
| 1 | 2019-01-01 02:00:00 | 13.80280 | 28.49950 | 110.0640 | 63.6111 | 55.9187 | 43.8401 | NaN | 33.6542 | 30.698 | ... | 73.4132 | 54.4664 | 19.0619 | 23.1494 | 10.7420 | 25.9358 | 47.9076 | NaN | NaN | 11.78830 |
| 2 | 2019-01-01 03:00:00 | 9.94056 | 11.12060 | 107.9410 | 48.3540 | 41.3488 | 22.8383 | NaN | 13.6030 | 28.262 | ... | 50.2355 | 50.4599 | 43.7717 | 21.0711 | 12.0391 | 24.5725 | 22.8309 | 5.57095 | NaN | 8.69917 |
| 3 | 2019-01-01 04:00:00 | 6.75889 | 5.57358 | 94.5489 | 34.6621 | 29.8697 | 20.1829 | NaN | 17.4302 | 26.522 | ... | 37.5872 | 34.8090 | 64.0139 | 21.1671 | 13.1849 | 20.6336 | 20.5900 | 5.77369 | NaN | 5.96861 |
| 4 | 2019-01-01 05:00:00 | 7.88722 | 6.56224 | 67.8800 | 14.2870 | 17.6000 | 18.7345 | NaN | 23.0878 | 24.260 | ... | 22.6446 | 30.6517 | 43.6111 | 21.0774 | 14.0005 | 19.4194 | 27.0838 | 6.15494 | NaN | 7.80778 |
5 rows × 64 columns
Zdecydowano się znaleźć dane dotyczące lokalizacji powyższych stacji pomiarowych, aby dodać uwzględnienie czynnika przestrzennego w dalszej analizie. Być może będą przydatne m.in. przy uzupełnianiu brakujących wartości lub przy sprawdzaniu zależności pomiędzy pomiarami w danym województwie. Zmieniono nazwy kolumn na takie, które odpowiadają nazwom stacji.
kody_stacji = pd.read_excel(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\Metadane oraz kody stacji i stanowisk pomiarowych (1).xlsx")
kody_stacji.head(5)
| Nr | Kod stacji | Kod międzynarodowy | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Data uruchomienia | Data zamknięcia | Typ stacji | Typ obszaru | Rodzaj stacji | Województwo | Miejscowość | Adres | WGS84 φ N | WGS84 λ E | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | DsBialka | NaN | Białka | NaN | 1990-01-03 | 2005-12-31 | przemysłowa | podmiejski | kontenerowa stacjonarna | DOLNOŚLĄSKIE | Białka | NaN | 51.197783 | 16.117390 |
| 1 | 2 | DsBielGrot | NaN | Bielawa - ul. Grota Roweckiego | NaN | 1994-01-02 | 2003-12-31 | tło | miejski | w budynku | DOLNOŚLĄSKIE | Bielawa | ul. Grota Roweckiego 6 | 50.682510 | 16.617348 |
| 2 | 3 | DsBogatFrancMOB | PL0602A | Bogatynia Mobil | DsBogatMob | 2015-01-01 | 2015-12-31 | tło | miejski | mobilna | DOLNOŚLĄSKIE | Bogatynia | ul. Francuska/Kręta | 50.940998 | 14.916790 |
| 3 | 4 | DsBogChop | PL0315A | Bogatynia - Chopina | NaN | 1996-01-01 | 2013-12-31 | przemysłowa | miejski | kontenerowa stacjonarna | DOLNOŚLĄSKIE | Bogatynia | ul. Chopina 35 | 50.905856 | 14.967175 |
| 4 | 5 | DsBogZatonieMob | PL0576A | Bogatynia - Mobil | NaN | 2012-01-01 | 2012-12-31 | przemysłowa | miejski | mobilna | DOLNOŚLĄSKIE | Bogatynia | ul. Konrada, Zatonie | 50.943245 | 14.913327 |
kody_stacji = kody_stacji[['Kod stacji','Nazwa stacji', 'Stary Kod stacji \n(o ile inny od aktualnego)', 'Typ obszaru', 'Województwo', 'Miejscowość', 'WGS84 φ N', 'WGS84 λ E']]
kody_stacji.rename(columns = {'WGS84 φ N': 'lat'}, inplace = True)
kody_stacji.rename(columns = {'WGS84 λ E': 'lon'}, inplace = True)
kody_stacji.head(5)
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 0 | DsBialka | Białka | NaN | podmiejski | DOLNOŚLĄSKIE | Białka | 51.197783 | 16.117390 |
| 1 | DsBielGrot | Bielawa - ul. Grota Roweckiego | NaN | miejski | DOLNOŚLĄSKIE | Bielawa | 50.682510 | 16.617348 |
| 2 | DsBogatFrancMOB | Bogatynia Mobil | DsBogatMob | miejski | DOLNOŚLĄSKIE | Bogatynia | 50.940998 | 14.916790 |
| 3 | DsBogChop | Bogatynia - Chopina | NaN | miejski | DOLNOŚLĄSKIE | Bogatynia | 50.905856 | 14.967175 |
| 4 | DsBogZatonieMob | Bogatynia - Mobil | NaN | miejski | DOLNOŚLĄSKIE | Bogatynia | 50.943245 | 14.913327 |
Wczytano dane znalezione na stronie GIOS, zawierające szczegółowe informacje o stacjach - ich lokalizację, województwo, typ obszaru itd. Wybrano kolumny, które uważano za istotne w dalszej analizie. Zmieniono nazwy kolumn na przyjemniejsze w użytku.
data_melted = data.melt(id_vars = 'Date', var_name = 'Kod stacji', value_name='PM2.5')
data_melted.head(10)
| Date | Kod stacji | PM2.5 | |
|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | DsDusznikMOB | 33.40530 |
| 1 | 2019-01-01 02:00:00 | DsDusznikMOB | 13.80280 |
| 2 | 2019-01-01 03:00:00 | DsDusznikMOB | 9.94056 |
| 3 | 2019-01-01 04:00:00 | DsDusznikMOB | 6.75889 |
| 4 | 2019-01-01 05:00:00 | DsDusznikMOB | 7.88722 |
| 5 | 2019-01-01 06:00:00 | DsDusznikMOB | 14.40720 |
| 6 | 2019-01-01 07:00:00 | DsDusznikMOB | 7.48833 |
| 7 | 2019-01-01 08:00:00 | DsDusznikMOB | 5.67083 |
| 8 | 2019-01-01 09:00:00 | DsDusznikMOB | 11.72390 |
| 9 | 2019-01-01 10:00:00 | DsDusznikMOB | 11.72530 |
Zdecydowano się przekształcić DataFrame 'data', tak aby usprawnić łączenie ze szczegółowymi kodami stacji. Teraz nazwy stacji nie są w osobnych kolumnach tylko w jednej połączonej.
data_merged = pd.merge(data_melted, kody_stacji, on = 'Kod stacji', how = 'left')
data_merged.head(5)
| Date | Kod stacji | PM2.5 | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | DsDusznikMOB | 33.40530 | Duszniki-Zdrój | NaN | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 1 | 2019-01-01 02:00:00 | DsDusznikMOB | 13.80280 | Duszniki-Zdrój | NaN | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 2 | 2019-01-01 03:00:00 | DsDusznikMOB | 9.94056 | Duszniki-Zdrój | NaN | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 3 | 2019-01-01 04:00:00 | DsDusznikMOB | 6.75889 | Duszniki-Zdrój | NaN | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 4 | 2019-01-01 05:00:00 | DsDusznikMOB | 7.88722 | Duszniki-Zdrój | NaN | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
Połączono dwa zestawy w jeden szczegółowy.
not_merged = data_merged[data_merged['Nazwa stacji'].isnull()]
not_merged['Kod stacji'].unique()
array(['LbNaleczow', 'MzKonJezMos', 'PdSuwPulaskp', 'PmGdaLeczk08',
'ZpSzczAndr01', 'ZpSzczPils02'], dtype=object)
6 nazw stacji nie znalzało swoich odpowiedników, więc zostaną podjęte próby ich odnalezienia, o ile istnieją
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'LbNaleczow']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
kody_stacji_notna = kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'].notna()]
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('LbNaleczow')]
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 296 | LbNaleczAlMa | Nałęczów, al. Małachowskiego | LbNaleczowMOB, LbNaleczow | podmiejski | LUBELSKIE | Nałęczów | 51.284931 | 22.210242 |
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'MzKonJezMos']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 536 | MzKonJezWieMOB | Konstancin-Jeziorna, ul. Wierzejewskiego | MzKonJezMos | podmiejski | MAZOWIECKIE | Konstancin-Jeziorna | 52.080625 | 21.111186 |
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'PdSuwPulaskp']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 688 | PdSuwPulask2 | Suwałki, ul. Pułaskiego 26 | PdSuwPulaskp | miejski | PODLASKIE | Suwałki | 54.115897 | 22.938464 |
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'PmGdaLeczk08']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('PmGdaLeczk08')]
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 799 | PmGdaLeczkow | Gdańsk, ul. Leczkowa | Pm.a08a, PmGdaLeczk08, PmGdaLecz08m | miejski | POMORSKIE | Gdańsk | 54.380279 | 18.620274 |
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'ZpSzczAndr01']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('ZpSzczAndr01')]
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 1082 | ZpSzczAndrze | Szczecin, ul. Andrzejewskiego | ZpSzczecin001, ZpSzczAndr01 | miejski | ZACHODNIOPOMORSKIE | Szczecin | 53.380975 | 14.663347 |
kody_stacji[kody_stacji['Stary Kod stacji \n(o ile inny od aktualnego)'] == 'ZpSzczPils02']
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
kody_stacji_notna[kody_stacji_notna['Stary Kod stacji \n(o ile inny od aktualnego)'].str.contains('ZpSzczPils02')]
| Kod stacji | Nazwa stacji | Stary Kod stacji \n(o ile inny od aktualnego) | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|
| 1090 | ZpSzczPilsud | Szczecin, ul. Piłsudskiego | ZpSzczecin002, ZpSzczPils02 | miejski | ZACHODNIOPOMORSKIE | Szczecin | 53.432169 | 14.5539 |
data_melted['Kod stacji'] = data_merged['Kod stacji'].replace({'LbNaleczow' : 'LbNaleczAlMa',
'MzKonJezMos': 'MzKonJezWieMOB',
'PdSuwPulaskp': 'PdSuwPulask2',
'PdSuwPulaskp': 'PdSuwPulask2',
'PmGdaLeczk08': 'PmGdaLeczkow',
'ZpSzczAndr01': 'ZpSzczAndrze',
'ZpSzczPils02': 'ZpSzczPilsud'})
Zatąpiono stare kody stacji nowymi, aby wszystkie dane się odpowiednio połączyły.
data_merged = pd.merge(data_melted, kody_stacji, on = 'Kod stacji', how = 'left')
data_merged = data_merged[['Date', 'Kod stacji', 'PM2.5', 'Nazwa stacji', 'Typ obszaru', 'Województwo', 'Miejscowość', 'lat', 'lon']]
data_merged.head(5)
| Date | Kod stacji | PM2.5 | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | DsDusznikMOB | 33.40530 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 1 | 2019-01-01 02:00:00 | DsDusznikMOB | 13.80280 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 2 | 2019-01-01 03:00:00 | DsDusznikMOB | 9.94056 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 3 | 2019-01-01 04:00:00 | DsDusznikMOB | 6.75889 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
| 4 | 2019-01-01 05:00:00 | DsDusznikMOB | 7.88722 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 |
not_merged = data_merged[data_merged['Nazwa stacji'].isnull()]
not_merged
| Date | Kod stacji | PM2.5 | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
Upewniono się, że wszystkie dane zostały połączone.
# poland = gpd.read_file(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\polska (1)\polska.shp")
woj = gpd.read_file(r"C:\Users\User\Desktop\magister\semestr1\ML\projekt_ocena_2\wojewodztwa\wojewodztwa.shp")
crs_polska = 'epsg:4326'
woj = woj.to_crs(crs_polska)
fig, ax = plt.subplots(figsize = (10, 10))
woj.plot(ax = ax, facecolor = 'lightblue', edgecolor = 'black')
ax.scatter(x = data_merged['lon'], y = data_merged['lat'], color = 'red', label = 'Stacje', s = 4)
plt.show()
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Kod stacji PmGdaLeczkow 53.059361 ZpSzczBudzWosMOB 18.493151 KpToruDziewu 17.659817 WmElbBazynsk 12.043379 KpBydPlPozna 9.691781 SkSkarZielnaMOB 9.143836 LdLodzGdansk 7.728311 PkHorZdrParkMOB 7.682648 LuZielKrotka 7.100457 WmGoldUzdrowMOB 6.929224 LuNowaSolMOB 6.803653 OpKKozBSmial 6.735160 SlKatoKossut 6.495434 LbNaleczAlMa 6.495434 MzWarWokalna 6.289954 PkPrzemGrunw 5.878995 LuWsKaziWiel 5.673516 LdZgieMielcz 5.616438 MzWarTolstoj 5.182648 KpMogiNowMOB 5.034247 Name: PM2.5, dtype: float64
Sprawdzono % brakujących wartości stężeń PM2.5 dla każdej ze stacji. Dla jednej stacji braki wynoszą ponad 50%. Braki w kolejnych stacjach są znacznie mniejsze. W 4 stacjach większe od 10%.
fig, axes = plt.subplots(1, 1, figsize = (20,10))
plt.plot(data_merged['Date'].unique(), data_merged['PM2.5'][data_merged['Kod stacji'] == 'PmGdaLeczkow'])
plt.xlabel('Date')
plt.ylabel('PM2.5')
plt.xlim(data_merged['Date'].min())
plt.title('PM2.5 values for station PmGdaLeczkow')
plt.show()
data_merged[data_merged['Województwo'] == 'PmGdaLeczkow']
| Date | Kod stacji | PM2.5 | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon |
|---|
Z racji braku innych obserwacji w województwie pomorskim i dużej ilości brakujących danych w początkowym okresie pomiarów zdecydowano się usunąć stację.
data_merged = data_merged[data_merged['Kod stacji'] != 'PmGdaLeczkow']
data_merged.groupby('Województwo')['Kod stacji'].nunique().sort_values()
Województwo WIELKOPOLSKIE 1 LUBELSKIE 2 OPOLSKIE 2 PODLASKIE 3 WARMIŃSKO-MAZURSKIE 3 ZACHODNIOPOMORSKIE 3 ŁÓDZKIE 3 ŚLĄSKIE 3 ŚWIĘTOKRZYSKIE 3 KUJAWSKO-POMORSKIE 4 LUBUSKIE 4 MAŁOPOLSKIE 4 DOLNOŚLĄSKIE 5 PODKARPACKIE 7 MAZOWIECKIE 15 Name: Kod stacji, dtype: int64
srednia_wojewodzka_godz = data_merged.groupby(['Date', 'Województwo'])['PM2.5'].transform('mean')
data_merged['PM2.5'] = data_merged['PM2.5'].fillna(srednia_wojewodzka_godz)
data_merged['srednia_wojewodzka_godz'] = srednia_wojewodzka_godz
data_merged['PM2.5'] = data_merged['PM2.5'].fillna(srednia_wojewodzka_godz)
data_merged.head(5)
| Date | Kod stacji | PM2.5 | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon | srednia_wojewodzka_godz | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | DsDusznikMOB | 33.40530 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 82.543420 |
| 1 | 2019-01-01 02:00:00 | DsDusznikMOB | 13.80280 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 54.379220 |
| 2 | 2019-01-01 03:00:00 | DsDusznikMOB | 9.94056 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 43.740992 |
| 3 | 2019-01-01 04:00:00 | DsDusznikMOB | 6.75889 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 34.282634 |
| 4 | 2019-01-01 05:00:00 | DsDusznikMOB | 7.88722 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 22.843292 |
Zdecydowano się zastąpić braki w danych średnią wojewódzką godzinową dla każdego z województw. Dane wciąż zawierają brakujące wartości ze względu na brak danych dla niektórych godzin.
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Kod stacji WpKaliSawick 2.420091 LdLodzGdansk 0.034247 LbLubObywate 0.034247 LbNaleczAlMa 0.034247 LdLodzCzerni 0.034247 LdZgieMielcz 0.034247 MzWarTolstoj 0.011416 MzPlocMiReja 0.011416 MzRadTochter 0.011416 MzSiedKonars 0.011416 MzWarAlNiepo 0.011416 SkKielTargow 0.011416 MzWarChrosci 0.011416 MzWarKondrat 0.011416 MzWarBajkowa 0.011416 MzWarWokalna 0.011416 MzOtwoBrzozo 0.011416 MzZyraRoosev 0.011416 OpKKozBSmial 0.011416 OpPrudPodgor 0.011416 Name: PM2.5, dtype: float64
srednia_dzienna = data_merged.groupby('Date')['PM2.5'].transform('mean')
maska_braki = (data_merged['Kod stacji'] == 'WpKaliSawick') & data_merged['PM2.5'].isna()
data_merged.loc[maska_braki, 'PM2.5']
indeksy = data_merged.loc[maska_braki, 'PM2.5'].index
# np.sum(srednia_dzienna_bez_WpKaliSawick.index == 516857)
srednia_dzienna.loc[indeksy]
data_merged.loc[maska_braki, 'PM2.5'] = srednia_dzienna.loc[indeksy]
Dla kodu stacji 'WpKaliSawick' wybrano metodę zastępowania braków średnią ze wszystkich województw.
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)[0:20]
Kod stacji LdLodzGdansk 0.034247 LbLubObywate 0.034247 LbNaleczAlMa 0.034247 LdLodzCzerni 0.034247 LdZgieMielcz 0.034247 MzWarTolstoj 0.011416 MzPiasPulask 0.011416 MzPlocMiReja 0.011416 MzRadTochter 0.011416 MzSiedKonars 0.011416 MzWarAlNiepo 0.011416 SkKielTargow 0.011416 MzWarChrosci 0.011416 MzWarKondrat 0.011416 MzWarBajkowa 0.011416 MzMinMazKaziMOB 0.011416 MzWarWokalna 0.011416 MzZyraRoosev 0.011416 OpKKozBSmial 0.011416 OpPrudPodgor 0.011416 Name: PM2.5, dtype: float64
x = data_merged[data_merged['Kod stacji'] == 'LdZgieMielcz']['PM2.5']
x[x.isna()]
117434 NaN 117435 NaN 121056 NaN Name: PM2.5, dtype: float64
data_merged['PM2.5'] = data_merged.groupby('Kod stacji')['PM2.5'].apply(lambda x: x.fillna(method='ffill', limit=2).fillna(method='bfill', limit=2))
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3740211283.py:1: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use
>>> .groupby(..., group_keys=False)
To adopt the future behavior and silence this warning, use
>>> .groupby(..., group_keys=True)
data_merged['PM2.5'] = data_merged.groupby('Kod stacji')['PM2.5'].apply(lambda x: x.fillna(method='ffill', limit=2).fillna(method='bfill', limit=2))
Dla powyższych stacji braki występują tylko dla jednej, bądź dla trzech godzin. Jeśli występują dla trzech godzin to sprawdzono, że co najwyżej dwa braki występują po sobie, więc zdecydowano się je zastąpić średnią z 2 godzin poprzedzających i następujących po braku.
(data_merged['PM2.5'].isna().groupby(data_merged['Kod stacji']).mean() * 100).sort_values(ascending=False)
Kod stacji
DsDusznikMOB 0.0
PkPrzemGrunw 0.0
MzWarKondrat 0.0
MzWarTolstoj 0.0
MzWarWokalna 0.0
...
MzOtwoBrzozo 0.0
MzPiasPulask 0.0
MzPlocMiReja 0.0
MzRadTochter 0.0
ZpSzczPilsud 0.0
Name: PM2.5, Length: 62, dtype: float64
Sprawdzenie,czy wszystkie dane nie zawierają już brakujących wartości PM2.5.
plt.figure(figsize=(12, 6))
for kod_stacji in data_merged['Kod stacji'].unique():
dane_stacji = data_merged[data_merged['Kod stacji'] == kod_stacji]
plt.plot(dane_stacji['Date'], dane_stacji['PM2.5'])
plt.title('PM2.5 dla każdej stacji')
plt.xlabel('Data i godzina')
plt.ylabel('PM2.5')
plt.xticks(rotation=45)
plt.show()
data_merged['Date'] = pd.to_datetime(data_merged['Date'])
data_merged['day'] = data_merged['Date'].dt.dayofweek
data_merged['month'] = data_merged['Date'].dt.month
data_merged['weekend'] = data_merged['day'].isin([5, 6]).astype(int)
data_merged['year'] = data_merged['Date'].dt.year
def get_season(month):
if month in [3, 4, 5]:
return 0
elif month in [6, 7, 8]:
return 1
elif month in [9, 10, 11]:
return 2
else:
return 3
data_merged['season'] = data_merged['month'].apply(get_season)
Zdecydowano się dodać nowe cechy takie jak: dzień, miesiąc, rok oraz informację czy jest weekend czy dzień roboczy.
unique_station_ids = data_merged['Kod stacji'].unique()
random_station_ids = random.sample(list(unique_station_ids), 1)
print(random_station_ids)
['LbLubObywate']
stacja = data_merged[(data_merged['month'] == 1) & (data_merged['year'] == 2019)]
stacja = stacja[stacja['Kod stacji'] == 'WmOlsPuszkin']
stacja = stacja.set_index("Date")
statistics = ["mean", "std", "max"]
windows = ['1h', '3h', '6h', '12h','1D', '3D']
for window in windows:
rolling_df = stacja['PM2.5'].rolling(window = window, closed='left').agg(statistics)
fig, axes = plt.subplots(len(statistics), 1, figsize = [20, 16])
fig.suptitle(f'Okno = {window}')
for idx, stat in enumerate(statistics):
sns.lineplot(data=stacja['PM2.5'], ax = axes[idx], color='red', label='PM2.5')
sns.lineplot(data=rolling_df[stat], ax=axes[idx], color='blue', label=f'{stat} PM2.5')
axes[idx].set_ylabel(f'{stat.capitalize()} PM2.5')
plt.show()
# stacja = data_merged[(data_merged['month'] == 1) & (data_merged['year'] == 2019)]
data_stacja = data_merged[data_merged['Kod stacji'] == 'WmOlsPuszkin']
data_stacja = data_stacja.set_index("Date")
statistics = ["mean", "std", "max"]
windows = ['7D', '30D']
for window in windows:
rolling_df = data_stacja['PM2.5'].rolling(window = window, closed='left').agg(statistics)
fig, axes = plt.subplots(len(statistics), 1, figsize = [20, 16])
fig.suptitle(f'Okno = {window}')
for idx, stat in enumerate(statistics):
sns.lineplot(data=data_stacja['PM2.5'], ax = axes[idx], color='red', label='PM2.5')
sns.lineplot(data=rolling_df[stat], ax=axes[idx], color='blue', label=f'{stat} PM2.5')
axes[idx].set_ylabel(f'{stat.capitalize()} PM2.5')
plt.show()
Od razu widać, że jednogodzinne okno nie niesie żadnej wartości, analogicznie dla 3godzinnego oraz 6godzinnego okna. Średnia dla 12 godzin ma już jakikolwiek sens. Zdecydowano się wybrać również okno ze średnią 1dniową oraz maximum 1dniowe. Wybrano również maxmium 30dniowe.
data_merged.set_index('Date', inplace=True)
def calculate_rolling_stats(group):
group['rolling_mean_12_h'] = group['PM2.5'].rolling(window='12h', closed='left').mean()
group['rolling_mean_1D'] = group['PM2.5'].rolling(window='1D', closed='left').mean()
group['rolling_max_1D'] = group['PM2.5'].rolling(window='1D', closed='left').max()
group['rolling_max_30D'] = group['PM2.5'].rolling(window='30D', closed='left').max()
return group
data_with_rolling_stats = data_merged.groupby('Kod stacji').apply(calculate_rolling_stats)
data_with_rolling_stats.reset_index(inplace=True)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\2148299505.py:10: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use
>>> .groupby(..., group_keys=False)
To adopt the future behavior and silence this warning, use
>>> .groupby(..., group_keys=True)
data_with_rolling_stats = data_merged.groupby('Kod stacji').apply(calculate_rolling_stats)
Dodano do zbioru danych cechy tj. średnia 12 godzinna, średnia 1 dniowa, maximum 1 dniowe oraz maximum 30 dniowe.
unique_station_ids = data_merged['Kod stacji'].unique()
random_station_ids = random.sample(list(unique_station_ids), 5)
print(random_station_ids)
['KpBydPlPozna', 'WmGoldUzdrowMOB', 'MzWarChrosci', 'PkMielBierna', 'MzZyraRoosev']
data_jan_2020 = data_merged[(data_merged['year'] == 2019) & (data_merged['month'] == 1)]
for i in random_station_ids:
pm25 = data_jan_2020[data_jan_2020['Kod stacji'] == i]['PM2.5']
plot_pacf(pm25, lags=24)
plt.show()
Dla większości przypadów najsilniejszą korelację osiągają opóźnienia 1godzinne, 2godzinne, 3godzinne, 20godzinne oraz 24 godzinne. Dlatego też zdecydowano się takie lagi dodać.
def calculate_rolling_stats(group):
lag1 = group['PM2.5'].shift(1)
group['lag1_h'] = lag1
lag2 = group['PM2.5'].shift(2)
group['lag2_h'] = lag2
lag3 = group['PM2.5'].shift(3)
group['lag3_h'] = lag3
lag20 = group['PM2.5'].shift(20)
group['lag20_h'] = lag20
lag24 = group['PM2.5'].shift(24)
group['lag24_h'] = lag24
return group
data_with_lags = data_with_rolling_stats.groupby('Kod stacji').apply(calculate_rolling_stats)
data_with_lags.reset_index(inplace=True)
data_with_lags = data_with_lags.drop(['index'], axis = 1)
C:\Users\User\AppData\Local\Temp\ipykernel_23260\2090405549.py:18: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use
>>> .groupby(..., group_keys=False)
To adopt the future behavior and silence this warning, use
>>> .groupby(..., group_keys=True)
data_with_lags = data_with_rolling_stats.groupby('Kod stacji').apply(calculate_rolling_stats)
data_with_lags.head(5)
| Date | Kod stacji | PM2.5 | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon | srednia_wojewodzka_godz | ... | season | rolling_mean_12_h | rolling_mean_1D | rolling_max_1D | rolling_max_30D | lag1_h | lag2_h | lag3_h | lag20_h | lag24_h | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2019-01-01 01:00:00 | DsDusznikMOB | 33.40530 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 82.543420 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 2019-01-01 02:00:00 | DsDusznikMOB | 13.80280 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 54.379220 | ... | 3 | 33.405300 | 33.405300 | 33.4053 | 33.4053 | 33.40530 | NaN | NaN | NaN | NaN |
| 2 | 2019-01-01 03:00:00 | DsDusznikMOB | 9.94056 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 43.740992 | ... | 3 | 23.604050 | 23.604050 | 33.4053 | 33.4053 | 13.80280 | 33.40530 | NaN | NaN | NaN |
| 3 | 2019-01-01 04:00:00 | DsDusznikMOB | 6.75889 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 34.282634 | ... | 3 | 19.049553 | 19.049553 | 33.4053 | 33.4053 | 9.94056 | 13.80280 | 33.4053 | NaN | NaN |
| 4 | 2019-01-01 05:00:00 | DsDusznikMOB | 7.88722 | Duszniki-Zdrój | miejski | DOLNOŚLĄSKIE | Duszniki-Zdrój | 50.402645 | 16.393319 | 22.843292 | ... | 3 | 15.976887 | 15.976887 | 33.4053 | 33.4053 | 6.75889 | 9.94056 | 13.8028 | NaN | NaN |
5 rows × 24 columns
data_with_lags.set_index('Date', inplace = True)
data_with_lags.sort_index(inplace = True)
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split, GridSearchCV
import xgboost as xgb
from sklearn.metrics import mean_squared_error, r2_score, mean_squared_error
from sklearn.model_selection import TimeSeriesSplit
X = data_with_lags.drop(['PM2.5'], axis=1)
y = data_with_lags['PM2.5']
columns_to_encode = [0, 1, 2, 3, 4]
for c in columns_to_encode:
labelencoder = LabelEncoder()
X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)` X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c]) C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)` X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c]) C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)` X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c]) C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)` X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c]) C:\Users\User\AppData\Local\Temp\ipykernel_23260\3553582975.py:5: DeprecationWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)` X.iloc[:, c] = labelencoder.fit_transform(X.iloc[:, c])
Zenkodowano kolumny kategoryczne, aby móc użyć je później w modelu.
X.head(5)
| Kod stacji | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon | srednia_wojewodzka_godz | day | month | ... | season | rolling_mean_12_h | rolling_mean_1D | rolling_max_1D | rolling_max_30D | lag1_h | lag2_h | lag3_h | lag20_h | lag24_h | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Date | |||||||||||||||||||||
| 2019-01-01 01:00:00 | 0 | 4 | 0 | 0 | 4 | 50.402645 | 16.393319 | 82.543420 | 1 | 1 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2019-01-01 01:00:00 | 5 | 3 | 0 | 1 | 3 | 53.121764 | 17.987906 | 46.479233 | 1 | 1 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2019-01-01 01:00:00 | 58 | 11 | 0 | 10 | 12 | 51.747950 | 18.049063 | 75.200000 | 1 | 1 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2019-01-01 01:00:00 | 43 | 7 | 1 | 7 | 8 | 50.192100 | 23.361425 | 40.202483 | 1 | 1 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2019-01-01 01:00:00 | 24 | 22 | 0 | 4 | 21 | 52.178766 | 21.560968 | 45.312354 | 1 | 1 | ... | 3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 22 columns
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle = False)
fig, axes = plt.subplots(1,1, figsize = (20,10))
plt.plot(y_train, color = 'red')
plt.plot(y_test, color = 'blue')
[<matplotlib.lines.Line2D at 0x14fafb94b90>]
Podzielono zbiór danych na treningowy i testowy w stosunku 75% do 25%. Zdecydowano się na wykonanie predykcji na danych ogólnych (nie dla każdej stacji z osobna) ze względu, że jest to bardziej optymalne obliczeniowo i mamy jeden model, który możemy zastosować dla całości danych, a nie musimy zależnie od stacji używać innego modelu.
model = xgb.XGBRegressor()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
rmse = mean_squared_error(y_test, y_pred, squared=False)
print("RMSE:", rmse)
RMSE: 6.247570667847238
Błąd średniokwadratowy dla modelu wynosi około 6.25% przy domyślnych parametrach.
PLUSY:
+zawiera techniki regularyzacji, które pomagają zapobiec nadmiernemu dopasownaniu,
+wysoka dokładność,
+dobrze radzi sobie z dużymi datasetami,
+nie ma potrzeby normalizowania danych, radzi sobie z brakującymi wartościami
MINUSY:
-wymaga dostosowania hiperparametróww, przekonano się, że jest to bardzo czasochłonne,
-pomimo technik regularyzacji algorytm nadal może przeuczyć dane treningowe,
-interpretacja i zrozumienie w jaki sposób algorytm osiąga przewidywania jest trudna,
-trudne do dostrojenia ze względu na dużą ilość hiperparametrów.
param_grid = {
'max_depth': [3, 5, 7],
'learning_rate': [0.1, 0.01, 0.001],
'n_estimators': [100, 200, 300],
}
grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)
print("Najlepsze parametry:", grid_search.best_params_)
best_model = grid_search.best_estimator_
y_pred = best_model.predict(X_test)
Najlepsze parametry: {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 100}
Po wykonaniu tuningu parametrów najlepszymi okazują się być:
-'learning_rate': 0.1,
-'max_depth': 5,
-'n_estimators': 100
Próbowano wykonywać kolejne tuningi z większą liczbą parametrów, lecz procesy te były zbyt czasochłonne.
Hiperparametry to zewnętrzne zmienne konfiguracyjne, których Data Scientiści używają do zarządzania szkoleniem modeli uczenia maszynowego. Są ustawiane ręcznie przed szkoleniem modelu. Hiperparametry określają kluczowe funkcje, takie jak architektura modelu, szybkość uczenia się i złożoność modelu. Mają znaczący wpływ na dokładność i skuteczność modelu.
Strojenie hiperparametrów to proces wybierania optymalnych wartości hiperparametrów modelu uczenia maszynowego. Celem dostrajania hiperparametrów jest znalezienie wartości, które prowadzą do najlepszej wydajności danego zadania.
best_params = {'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 100}
final_model = xgb.XGBRegressor(**best_params)
final_model.fit(X_train, y_train)
y_pred = final_model.predict(X_test)
rmse = mean_squared_error(y_test, y_pred, squared=False)
print("RMSE:", rmse)
RMSE: 6.198891964836731
Po dostrojeniu hiperparametrów błąd średniokwadratowy zmalał w porównaniu do modelu z wartościami domyślnymi.
testing_data = pd.DataFrame(data=X_test, columns=X.columns)
testing_data['PM2.5'] = y_test
testing_data['pm25_przewidywane'] = y_pred
testing_data.head()
| Kod stacji | Nazwa stacji | Typ obszaru | Województwo | Miejscowość | lat | lon | srednia_wojewodzka_godz | day | month | ... | rolling_mean_1D | rolling_max_1D | rolling_max_30D | lag1_h | lag2_h | lag3_h | lag20_h | lag24_h | PM2.5 | pm25_przewidywane | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Date | |||||||||||||||||||||
| 2019-10-01 19:00:00 | 33 | 48 | 0 | 4 | 41 | 52.290864 | 21.042458 | 7.482966 | 1 | 10 | ... | 4.482526 | 6.322227 | 51.180685 | 5.011755 | 4.507727 | 4.507727 | 4.709338 | 1.987589 | 5.213366 | 5.518918 |
| 2019-10-01 19:00:00 | 26 | 28 | 1 | 4 | 27 | 52.191728 | 20.837489 | 7.482966 | 1 | 10 | ... | 5.584167 | 7.210000 | 260.490000 | 6.190000 | 5.400000 | 5.510000 | 5.060000 | 4.610000 | 7.100000 | 6.485516 |
| 2019-10-01 19:00:00 | 58 | 11 | 0 | 10 | 12 | 51.747950 | 18.049063 | NaN | 1 | 10 | ... | 6.613928 | 11.300000 | 65.493700 | 4.094110 | 4.182160 | 2.128030 | 8.866600 | 5.002220 | 7.516842 | 7.554872 |
| 2019-10-01 19:00:00 | 19 | 16 | 0 | 5 | 16 | 50.010575 | 19.949189 | 9.795942 | 1 | 10 | ... | 5.686784 | 15.587000 | 52.089600 | 6.675240 | 3.294000 | 3.873900 | 4.179670 | 8.146100 | 5.802210 | 7.633582 |
| 2019-10-01 19:00:00 | 9 | 20 | 0 | 2 | 19 | 51.259431 | 22.569133 | 8.867390 | 1 | 10 | ... | 4.987500 | 8.800000 | 64.600000 | 8.800000 | 7.700000 | 5.600000 | 4.400000 | 2.900000 | 11.900000 | 9.025276 |
5 rows × 24 columns
testing_data_stacja = testing_data[testing_data['Kod stacji'] == 33]
r2 = r2_score(y_test, y_pred)
print(f'R-squared: {r2}')
plt.figure(figsize=(12, 8))
plt.plot(testing_data_stacja['PM2.5'], label='Rzeczywiste PM2.5')
plt.plot(testing_data_stacja['pm25_przewidywane'], label='Przewidywane PM2.5')
plt.legend()
plt.title('Rzeczywiste vs przewidywane wartości PM2.5')
plt.xlabel('Index')
plt.ylabel("Stężenie PM2.5 (µg/m³)")
plt.show()
R-squared: 0.8744843687608252
# param_grid = {
# 'max_depth': [3, 5, 7],
# 'learning_rate': [0.05, 0.1, 0.2],
# 'n_estimators': [100, 200, 300],
# 'colsample_bytree': [0.6, 0.8, 1.0],
# 'subsample': [0.6, 0.8, 1.0],
# 'gamma': [0, 0.1, 0.2],
# 'reg_alpha': [0, 0.1, 0.2],
# 'reg_lambda': [0, 0.1, 0.2],
# 'min_child_weight': [1, 3, 5],
# 'max_delta_step': [0, 1, 2]
# }
# grid_search = GridSearchCV(estimator=xgb.XGBRegressor(), param_grid=param_grid, cv=3, scoring='neg_mean_squared_error')
# grid_search.fit(X_train, y_train)
# print("Najlepsze parametry:", grid_search.best_params_)
# best_model = grid_search.best_estimator_
# y_pred = best_model.predict(X_test)
# rmse = mean_squared_error(y_test, y_pred, squared=False)
# print("RMSE:", rmse)
tscv = TimeSeriesSplit(n_splits=5)
Walidacja Walk Forward (WFV) to technika walidacji krzyżowej szeregów czasowych stosowana do oceny wydajności modeli predykcyjnych. Walk forward polega na podziale zbioru danych na sekwencje ciągłych okien czasowych. Jest to szczególnie przydatne w przypadku danych uporządkowanych w czasie, gdzie sekwencja danych ma znaczenie. WFV zaprojektowano tak, aby był bardziej realistyczny w ocenie, jak dobrze model będzie generalizował na przyszłe, niewidoczne dane.
DZIAŁANIE:
Początkowy okres szkolenia i testu.
Trenowanie modelu.
Model testowy.
Przesunięcie okna.
Powtórzenie.
rmse_scores = []
for train_index, test_index in tscv.split(X_train):
X_train_fold, X_test_fold = X_train.iloc[train_index], X_train.iloc[test_index]
y_train_fold, y_test_fold = y_train.iloc[train_index], y_train.iloc[test_index]
print(f"Długości danych treningowych: {len(y_train_fold)}")
model.fit(X_train_fold, y_train_fold)
y_pred = model.predict(X_test_fold)
rmse = mean_squared_error(y_test_fold, y_pred, squared=False)
rmse_scores.append(rmse)
mean_rmse = sum(rmse_scores) / len(rmse_scores)
print("Mean RMSE:", mean_rmse)
Długości danych treningowych: 67890 Długości danych treningowych: 135780 Długości danych treningowych: 203670 Długości danych treningowych: 271560 Długości danych treningowych: 339450 Mean RMSE: 3.8984356076768236
Backtesting polega na ocenie wydajności modelu na danych historycznych, które nie zostały użyte do trenowania, symulując prawdziwe warunki na rynku lub środowisku. Celem backtestingu jest sprawdzenie, jak dobrze model radzi sobie ze zrozumieniem i przewidywaniem rzeczywistych zdarzeń na podstawie historycznych danych.
Testowanie wydajności modelu dla danych niebędących szeregami czasowymi odbywa się na zbiorze danych, który jest jednorodny w czasie, ale nie ma konkretnego związku z określonymi punktami w czasie. Celem testowania wydajności modelu dla danych niebędących szeregami czasowymi jest ocena ogólnej zdolności modelu do generalizacji na nowych, nieznanych danych, które nie mają struktury czasowej.
Podsumowując, backtesting ma specyficzny charakter związany z oceną modeli predykcyjnych w kontekście historycznych danych czasowych, podczas gdy testowanie wydajności modelu dla danych niebędących szeregami czasowymi koncentruje się na ocenie zdolności modelu do generalizacji na różne rodzaje danych bez struktury czasowej.
feature_importance = final_model.feature_importances_
feature_importance_df = pd.DataFrame({'Feature': X_train.columns, 'Importance': feature_importance})
feature_importance_df = feature_importance_df.sort_values(by='Importance', ascending=False)
plt.figure(figsize=(10, 6))
plt.barh(feature_importance_df['Feature'], feature_importance_df['Importance'])
plt.xlabel('Importance')
plt.ylabel('Feature')
plt.title('Feature Importance')
plt.show()
feature_importance_df
| Feature | Importance | |
|---|---|---|
| 17 | lag1_h | 0.836431 |
| 7 | srednia_wojewodzka_godz | 0.048263 |
| 15 | rolling_max_1D | 0.030158 |
| 16 | rolling_max_30D | 0.009065 |
| 1 | Nazwa stacji | 0.008964 |
| 19 | lag3_h | 0.007767 |
| 8 | day | 0.007317 |
| 18 | lag2_h | 0.005483 |
| 4 | Miejscowość | 0.005044 |
| 6 | lon | 0.004881 |
| 5 | lat | 0.004796 |
| 21 | lag24_h | 0.004277 |
| 20 | lag20_h | 0.004110 |
| 12 | season | 0.003913 |
| 13 | rolling_mean_12_h | 0.003774 |
| 14 | rolling_mean_1D | 0.003760 |
| 2 | Typ obszaru | 0.003722 |
| 0 | Kod stacji | 0.003030 |
| 9 | month | 0.002894 |
| 3 | Województwo | 0.002350 |
| 10 | weekend | 0.000000 |
| 11 | year | 0.000000 |
Można zaobserować (tak jak można było się spodziewać), że największą istotność ma kolumna z wartościami jednogodzinnego opóźnienia. Kolejnymi czynnikami, które mają jakiekolwiek znaczenie są cechy stworzone przy okazji podziału na okna czasowe oraz w trakcie tworzenia opóźnień.